program VMD

  use fft_wrap, only   : fft_init
  use mpi_wrap
  use efield
  use evolve

  use dimensions !, only : jx, jz, jv, jDir

  use parameters, only : restart, nt, dt, iDiag, iData, vPhase, runname
  use parameters, only : iTransAdvect, iEraseNonlin, iLinVlasov, iLinHydro, iFourier

  implicit none

  !-- major arrays ---------------------

  complex, allocatable  :: dfu(:, :, :, :, :)   ! distribution function
  complex, allocatable  :: uxMat(:, :, :, :)    ! velocity, x-component
  complex, allocatable  :: uyMat(:, :, :, :)    ! velocity, y-component, unused
  complex, allocatable  :: uzMat(:, :, :, :)    ! velocity, z-component, unused
  complex, allocatable  :: Ex2D(:, :, :)        ! electric field, x-component 
  complex, allocatable  :: Ey2D(:, :, :)        ! electric field, y-component 
  complex, allocatable  :: Ez2D(:, :, :)        ! electric field, z-component

  real,    allocatable  :: filter(:,:,:)        ! conical filter
 
  real,    allocatable  :: f0vec(:)             ! background distribution, gaussian
  complex, allocatable  :: fInit(:,:,:)         ! initial df  (1D equilibrium)


  !-- variables  ----------------------- 

  real       :: t                               ! time
  real       :: courant
  real       :: Courant_in_x                    ! advection Courant number 
  real       :: Courant_in_y                    ! advection Courant number 
  real       :: Courant_in_z                    ! advection Courant number 

  integer    :: it                              ! iteration count
  integer    :: iargc                           ! argument count

 
!-------------------------------------------------------------------------------------

  allocate ( dfu(jv, jz, jy, jx, jDir) )
  allocate ( uxMat(jz, jy, jx, jDir) )
  allocate ( uyMat(jz, jy, jx, jDir) )
  allocate ( uzMat(jz, jy, jx, jDir) )
  allocate ( Ex2D(jz, jy, jx) )
  allocate ( Ey2D(jz, jy, jx) )
  allocate ( Ez2D(jz, jy, jx) )

  allocate ( filter(jz,jy,jx) )

  allocate ( f0vec(jv) )
  allocate ( fInit(jv,jz,jDir) )


      !--  reading the run_name from the command line

      if(iargc().eq.0) then
         call getarg(0,runname)
         write(*,*)      'Usage: ',trim(runname),' <run name>'
         call mpi_abort(MPI_COMM_WORLD, -1, mpi_err)
         stop
      end if
      call getarg(1,runname)


      !-- set parameters, coordinate space, and the filter --

      call mpi_wrap_init()

      call read_input()

      call fft_init(dfu)

      call set_space_vectors()
      call set_filter(filter)

      call evolve_init()
      
      call efield_init()



      !-- restore from restart file or create initial conditions --

      call set_ic( dfu, uxMat, uyMat, f0vec, fInit, it)

      t  = it*dt

      call efield_evolve_E0(t)

      call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)


      !-- write initial diagnostics and binaries --

      if (.not. restart) then

         Courant_in_z = vPhase*dt/dz

         call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, t)

         call write_dispersion()
         call write_f0_initial(f0vec, dfu)

         call evolve_Uz( dfu, f0vec, uzMat)
         call write_bin(dfu, uxMat, uyMat, it, it/iData)

         write(msg, "('Output ' i4.4)") it/iData
         call mpi_wrap_msg

      end if

      write(msg,*) char(10), ' Begin temporal evolution', char(10)
      call mpi_wrap_msg


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------

 
do while (it .le. nt) 

   call dealias_5D(dfu)
   call dealias_4D(uxMat)
   call dealias_4D(uyMat)

   call evolve_F_long_advect(dfu, .false.)

   if (iTransAdvect) then
      call evolve_F_trans_advect( dfu, uxMat, uyMat, f0vec)
   end if


   call efield_evolve_E0(t)
   call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)

   call evolve_F_accel(  dfu, Ez2D, f0vec)
   call evolve_Uz(       dfu, f0vec, uzMat) 

   if (iTransAdvect) then
      call evolve_U_long_advect(  uxMat, uyMat,  uzMat)
      call evolve_U_trans_advect( uxMat, uyMat)
      call evolve_U_accel(        uxMat, uyMat, Ex2D, Ey2D)
      call filter_U_modulation(   uxMat, uyMat, filter)
   end if


   if (iEraseNonlin) then
      call evolve_F_erase_nonlin(dfu, fInit) 
      call evolve_U_erase_nonlin(uxMat, uyMat)
   end if

   t  = t + dt
   it = it + 1
   
   courant = maxval(abs(uzMat))*dt/dz
   call mpi_allreduce(courant, Courant_in_z, 1, &
        MPI_REAL8, MPI_MAX, MPI_COMM_WORLD, mpi_err)

   courant = maxval(abs(uxMat))*dt/dx
   call mpi_allreduce(courant, Courant_in_x, 1, &
        MPI_REAL8, MPI_MAX, MPI_COMM_WORLD, mpi_err)

   courant = maxval(abs(uyMat))*dt/dy
   call mpi_allreduce(courant, Courant_in_y, 1, &
        MPI_REAL8, MPI_MAX, MPI_COMM_WORLD, mpi_err)


   !-- diagnostics --

   if ( (it/iDiag)*iDiag.eq.it ) then      
   
     call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, t)

     !write(msg, "(' Diagnostcs, t =' e12.6)") t
     !call mpi_wrap_msg

   end if

   !-- output binary data for restart and post-processing --

   if ( (it/iData)*iData.eq.it ) then
     call write_bin(dfu, uxMat, uyMat, it, it/iData)
     write(msg, "(' Output ' i4.4)") it/iData
     call mpi_wrap_msg

   end if

   !-- check Courant conditions --
      
   if(  (Courant_in_z > 1.0) .and. .not. (iFourier .eq. 1) ) then
       write(msg,*) 'Courant for z advection =', Courant_in_z, '  is > 1, stop evolution'
       call mpi_wrap_msg
       go to 314
   endif
             
   if(  (Courant_in_x > 1.0)  .and. .not. (iLinVlasov .and. iLinHydro) ) then
       write(msg,*)   'Courant for x advection =', Courant_in_x, '  is > 1, stop evolution'
       call mpi_wrap_msg
       go to 314
   endif

   if(  (Courant_in_y > 1.0) .and. .not. (iLinVlasov .and. iLinHydro) ) then
      write(msg,*)   'Courant for y advection =', Courant_in_y, '  is > 1, stop evolution'
      call mpi_wrap_msg
      go to 314
   endif

end do      ! END TEMPORAL LOOP
 
314 continue



    write(msg,*) char(10), &
         ' Dvvvv should be > E_max / p_max^3 = ', maxval(Ez2D) / kVvec(jv/2)**3
    call mpi_wrap_msg

    write(msg,*) char(10), ' End temporal evolution', char(10)
    call mpi_wrap_msg

    !-- write some specrtral spectrum quantities --
       
    call mpi_wrap_exit()
    
    stop;

end

!-------------------------------------------------------------------------------------
